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The Landau paradigm of classifying phases by 
broken symmetries was demonstrated to be in- 
complete when it was realized that different quan- 
tum Hall states could only be distinguished by 
more subtle, topological properties [1]. Today, 
the role of topology as an underlying description 
of order has branched out to include topological 
band insulators, and certain featureless gapped 
Mott insulators with a topological degeneracy in 
the groundstate wavefunction. Despite intense fo- 
cus, very few candidates for these topologically 
ordered "spin liquids" exist. 

The main difficulty in finding systems that har- 
bour spin liquid states is the very fact that they 
violate the Landau paradigm, making conven- 
tional order parameters non-existent. Here, we 
uncover a spin liquid phase in a Bose-Hubbard 
model on the kagome lattice, and measure its 
topological order directly via the topological en- 
tanglement entropy. This is the first smoking-gun 
demonstration of a non-trivial spin liquid, identi- 
fied through its entanglement entropy as a gapped 
groundstate with emergent Z2 gauge symmetry. 

Quantum spin liquid phases are notoriously elusive, 
both in experimental materials and in theoretical mod- 
els. In part, this is due to the delicate balance of mi- 
croscopic interaction that must occur so that conven- 
tional symmetry-broken order is suppressed to low tem- 
peratures. The search is also hampered by the lack of 
a measurable order parameter, such as magnetization, 
that would offer a positive indicator of spin liquid be- 
havior. Instead, the current procedure of identifying 
spin liquids involves eliminating all possible order pa- 
rameters through exhaustive searches of correlation func- 
tions |2]. Theoretical work has nonetheless developed a 
classification scheme of gapped spin liquid states based 
on the topological degeneracy of their wavefunction. In 
fact, there exist several model Hamiltonians that have 
been proposed to contain spin liquid groundstates with 
the most trivial Z2 topological order (corresponding to 
a four- fold degeneracy on a torus). One of the earli- 
est was the triangular-lattice quantum dimer model [s) 
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FIG. 1: An L = 8 toroid (a periodic lattice with 8x8x3 
spins). In Equation O the subregions are Ai = aU6UcUd, 
A2 = a U 6 U c, A3 = a U 5 U (i, and = a U 6. The width of 
the annulus is = 6, and its thickness is r = 2. 



- where dimers are intended to be an effective descrip- 
tion of local singlet correlations, not physical spins. An- 
other paradigm in studying topological order has recently 
emerged in the toric code [4], as it is the simplest of a 
class of exactly solvable models (such as the Levin- Wen 
models [5]) describing different topological quantum field 
theories. Unfortunately, these models require somewhat 
artificial multi-spin interaction terms. 

In order for a model to be relevant for real physical sys- 
tems, it is essential to find a spin Hamiltonian with simple 
two-body interaction terms and a topologically ordered 
spin liquid state. One important step in this direction oc- 
curred with the introduction of a kagome Bose-Hubbard 
model by Balents, Girvin and Fisher [6]. Although this 
model has a four-spin interaction, it contains a Z2 spin 
liquid over an extended region of its phase diagram. Pre- 
viously, a related model was shown through quantum 
Monte Carlo (QMC) studies to have a featureless Mott 
insulating state f^J^^ where, like most experimental can- 
didates, the absence of order in correlation functions was 
the main evidence for spin liquid behavior. 

In 2006, Levin and Wen [9] and Kitaev and Preskill [lo| 
identified a quantity 27 called the topological entangle- 
ment entropy (EE) which is designed to replace the con- 
cept of an "order parameter" in a topologically ordered 
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system. Based on the idea that the spin liquid state is 
a type of collective paramagnet, the topological EE is 
designed to pick up non-local correlations in the ground- 
state wavefunction that are not manifest as conventional 
long-range order. However, these correlations contribute 
to the total entanglement between different subregions 
of the system A and its complement B (where AU B is 
the entire system). The EE between A and B can be 
quantified by the Renyi entropies, 

5„(A) = -^ln[Tr(pl)], (1) 

1 — n 

where pA is the reduced density matrix of region A. 
In a topologically ordered state, the non-local entangle- 
ment gives a topology-dependent subleading correction 
to "area-law" scaling of the EE of subregion A. In 2D, 

Sn{A) = ai — jj-l where a is a non-universal constant, 

£ is the boundary length between A and 5, and j is the 
number of disconnected boundary curves. In Levin and 
Wen's [9| construction (used in the calculations in this 
paper), the topological contribution can be isolated from 
the area-law scaling (plus any corner contributions) by 
considering separately the Renyi entropies on four 
differently shaped subregions (Figure [1]), 

27= lim [-5„(Ai) + 5„(A2) + 5„(A3)-^„(^4)]. 

(2) 

Naively, since calculating 7 requires complete knowl- 
edge of the groundstate wavefunction (through pa)^ pre- 
vious efforts to calculate it have been restricted to models 
which can be solved exactly either analytically (e.g. the 
toric code) or through numerical exact diagonalization 
on small size systems (e.g. the triangular lattice dimer 
model 0]). The ability to use 7 as a general tool to 
search for and characterize non-trivial topologically or- 
dered phases has been hindered by the inability to ac- 
cess the wavefunction in large-scale numerical methods, 
namely QMC, currently the only scalable quantum sim- 
ulation method in 2D and higher. However, with the 
recent introduction of measurement methods based on 
the replica tricky QMC is now able to access Sn{A) for 
n > 2 [13], therefore giving one a method to calculate 7 
in large-scale simulations of quantum spin liquids. 

Using Stochastic Series Expansion QMC [iBj, we 
simulate a hard-core Bose-Hubbard model on the kagome 
lattice, with nearest-neighbor hopping and a six-site po- 
tential around each lattice hexagon, 

H = -t^^lblbj + b.b]] + F^(no)^ (3) 

(ij) O 

where bj (bi) is the boson creation (annihilation) opera- 
tor, and no = X^^^Q^i, where n^ = blbi is the num- 
ber operator. As mentioned above, variations of this 
model with more complicated spin interactions are known 
to harbour a robust spin hquid groundstate In 



this paper, we consider the simplified Hamiltonian (j3j), 
with only nearest-neighbor hopping, which may be more 
amenable to construction for example in real cold atomic 
systems. We observe a transition at low temperature 
between a superfluid phase and an insulating phase for 
{V/t)c - 7.0665(15) (Figure [2]). For V/t > {V/T)c 
the superfluid density scales to zero, and density and 
bond correlators are featureless (similar to Ref. 7). This 
strongly suggests that the insulating phase is a spin liq- 
uid. To characterize it, we calculate the topological EE, 
Equation (|2]) with n = 2, which for a Z2 topological 
phase should approach 21n(2) in the limit T ^ 0. 
The regions Ai are shown in Figure [1] for an L = 8 sys- 
tem; these are scaled proportionally for the other systems 
sizes studied in this paper, where L is always a multiple 
of 8. Results for 7 as a function of inverse temperature 
/3 = t/T are shown in Figure [3] for several V/t. 

In the topological phase (V/t = 8) we see two dis- 
tinct plateaus, at differing temperatures, with a non-zero 
topological EE as T ^ 0. The phenomenon is known to 
occur in other models such as the toric code ^6j, where 
the topological EE at zero temperature of 2 ln(2) can be 
viewed as a sum of electric and magnetic contributions, 
each contributing ln(2). If the electric and magnetic de- 
fects have different energies, theory predicts two distinct 
plateaus corresponding to these individual crossover tem- 
peratures [l6|, as seen in our data. However, at any fixed 
non-zero temperature, in the limit of large L, the topolog- 
ical EE vanishes, as the probability of having thermally 
excited defects in the annulus Ai (Figure [1]) tends to 
unity. Indeed, under the assumption that the probabil- 
ity of having a defect is proportional to exp(— £^//cbT), 
where E is the defect energy, the temperature required 
to see accurate plateaus in the topological EE scales log- 
arithmically with L. In Figure O we show finite-size scal- 
ing data consistent with this logarithmic scaling. Note 
that our value for the topological EE at the higher-T 
plateau is indeed very close to ln(2), and becomes more 
accurately quantized at larger system sizes. The value of 
2 ln(2) at the lower-T plateau for L = 8 is not as accu- 
rately quantized, but is still approached. 

In the superfluid phase, the topological EE tends 
to zero as T ^ (Figure [3|). However, surprisingly, 
for V/t = 6 we observe a plateau in the topological 
EE at intermediate temperatures, T ^ t. One possi- 
ble explanation for this plateau can be understood by 
thinking of a simpler phase transition present in the 
toric code, induced by adding a parallel magnetic field. 
Consider a square-lattice toric code Hamiltonian H = 

-U E+ n,6+ Si - 9 En n.en Sf - h Ei Sf, where the 
first vertex term term penalizes vertices that do not have 
an even number of up spins on the legs of the neighboring 
bonds, and the second sum is over plaquettes. Suppose 
U ^ g. By increasing h/g^ we induce a phase transition 
from a topological phase to a trivial phase. In a non-zero 
temperature regime where U ^ T ^ g^ the problem be- 
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FIG. 2: (a) The schematic phase diagram of the kagome Bose-Hubbard model (|3)). Insets are the structure factor (Fourier 
transform of the equal-time density-density correlation function) for V/t — Q and 8 at high T (/3 = 3, top) and low T {V/t — 8 
and (5 — 48, bottom) in the spin liquid phase, (b) Data collapse of the superfluid density, which in the vicinity of a continuous 
phase transition should scale as ps(L, F/t, /3) = L~^F{\y/t — {V /t)c]L^^^ ^ ^ / L^). Here, F is the scaling function, z is the 
dynamical critical exponent, and v is the correlation length exponent. It follows from the above equation that if we plot psL 
as a function of \y/t — {V /t)c]L^^^ at fixed then the curves for different system sizes should collapse onto the universal 

curve F for appropriate values of v and (y/t)c, as shown for v — 0.6717, {y/t)c — 7.0665, and P/L — 2. 



comes classical: the quantum dynamics induced by g can 
be ignored and the vertex term restricts us to states de- 
scribed by closed loops of up spins. In this classical prob- 
lem, there is a phase transition as /i/T increases from a 
topological phase with long loops to a trivial phase with 
only short loops (this transition is dual to the 2D Ising 
transition). Thus, at high temperatures (small h/T) we 
see a topological EE, while at lower T the topological 
EE disappears. In fact, this kind of physics has been 
suggested to occur as a "cooperative paramagnet" in a 
related kagome lattice model ||8|]. 

The QMC results indicate a quantum critical point, 
separating the superfluid and topological phases, for a 
critical V/t located precisely by studying the finite size 
scaling of the superfluid density ps (Figure [2]). The data 
scales very well with the dynamical exponent z = 1 and 
the XY value [l^ for the correlation length exponent 
V = 0.6717 at {V/t)c = 7.0665(15). At T = near the 
quantum critical point, we expect good quantization of 
the topological EE whenever L is sufficiently large com- 
pared to the correlation length ^, so the topological EE 
may be controlled by a scaling function of L/^. At T > 0, 
the scaling of the topological EE in the quantum criti- 
cal fan appears not to have been considered previously; 
the plateau at intermediate temperatures for V/t = 6 
can perhaps also be understood as a manifestation of in- 
creasing T moving one from the zero temperature trivial 
phase into the quantum critical fan. Even at T = 0, the 
behavior of constant terms in the entropy at a critical 
point is largely unexplored, and may depend sensitively 
upon the geometry used to define it [22]. Scaling pre- 
dicts that near the critical point the topological EE is 
a function of and L/^, implying that in the topo- 



logical phase, corrections to 2 ln(2) should depend upon 
(L/^)^ exp(— /3/^), consistent with a defect energy of or- 
der 1/^. Based on the intermediate temperature plateau 
at F/t = 6, it seems likely that only one type of defect, 
the magnetic defect, becomes gapless at criticality. 

Discussion - In order to identify a topological phase, 
it is essential to perform non-local probes. Experimen- 
tally, such non-local probes could involve braiding oper- 
ations as in the proposal of In this paper, we have 
shown that the topological entanglement entropy (EE), 
calculated by QMC using the replica trick, is a practical 
numerical non-local probe. Other probes might be pos- 
sible, such as calculating the ground state degeneracy on 
lattices of different topology. However, that probe suffers 
from two drawbacks. First, studying surfaces of different 
Euler characteristic requires introducing defects, which 
is undesirable [23]. Second, while it would be possible in 
QMC to calculate the ground state degeneracy by inte- 
grating the specific heat, it requires accurate simulations 
at a temperature low enough to suppress all excitations 
- a regime where simulation ergodicity typically becomes 
a problem. In contrast, we have demonstrated that mea- 
surement of topological EE yields accurate quantization 
once one has suppressed all topologically non-trivial de- 
fects - and indeed one even sees accurate quantization 
at higher temperature where only one kind of defect is 
suppressed. Thus, we expect that replica QMC measure- 
ments of topological EE will be a fundamental technique 
in the characterization of non-trivial topological phases 
in the future. 
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FIG. 3: (a) The topological entanglement entropy (EE) measured on an L = 8 system as a function of inverse temperature 
(3 = t/T. The plateaus are a measure of the total quantum dimension 9], and should be ln(2) and 21n(2) for a Z2 spin liquid, 
(b) The approach of the topological EE to the first plateau, for different system sizes. The value of the crossover temperature 
(/3x), measured at 27 = ln(2)/2, shows a logarithmic dependence on system size (inset). 



HP^C initiative. Simulations were performed on the Bru- 
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